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Abstract 

We show that symplectic Runge-Kutta methods provide effective sym- 
plectic integrators for Hamiltonian systems with index one constraints. 
These include the Hamiltonian description of variational problems sub- 
ject to position and velocity constraints nondegenerate in the velocities, 
such as those arising in subRiemannian geometry. 



1 Introduction: constrained Hamiltonian systems 

We consider constrained Hamiltonian systems of the form 

Jz = VH{z), zeCcl" 1 (1) 

where z G M m , u:=\dzhJdz is a closed 2-forrrH H : W n -> R is a Hamiltonian, 
and C is a constraint submanifold such that i*oj (where i: C —> M. m is the 
inclusion of C in ]R m ) is nondegenerate, i.e., such that (C,i*u>) is a symplectic 
manifold. The dynamics on C depends only on the restricted Hamiltonian 
i* H and restricted symplectic form i*uj. Systems with holonomic (position) 
constraints take this form, with z = (q,p), lo = dq A dp, and C = {{q,p) : 
hi(q) = 0, Dhi(q)H p (q,p) = 0, 1 < i < k} consisting of primary and secondary 
constraints; a nondegeneracy assumption ensures that C is symplectic. The 
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We use vector notation in wedge products, writing dqAdp for ^ZI^Li dqiAdpi and dzhJdz 
for X]7j=i Jijd-Zi A dzj, where the dimension m is determined from the context. 
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widely used rattle method [31 [5] provides a (class of) symplectic integrators 
for this case when J is constant: it integrates in coordinates z with Lagrange 
multipliers to enforce the constraints. However, there are no known symplectic 
integrators for general constrained Hamiltonian systems of the form of Eq. ([!} . 

In this paper we describe a class of symplectic integrators for a class of 
Hamiltonian systems of the form ([!} containing constraints that can depend 
on both position and velocitjj^) The class includes the Hamiltonian description 
of problems arising from variational problems in subRiemannian geometry, in 
which velocities are constrained to lie in a given (nonintegrable) distribution. We 
give this application first. In the following proposition, the linear independence 
assumption on the constraints is equivalent to constraining the velocities to lie 
in a fc-dimensional distribution of the tangent space of the positions. 

Proposition 1. Let M be a symmetric nonsingular nxn mass matrix, V : K n — ► 
K a smooth potential, gi : R™ — > W l , i = l,...,k be smooth functions whose 
values are linearly independent for all arguments, and q be a smooth extremal 
with fixed endpoints for the functional 

S(q) = jf 1 L(t, q, q)dt = jf ' (^q T Mq - V(q)^j dt (2) 

subject to the constraints gi(q) ■ q = 0, i = 1, . . . , k. Then 

Jz = VH{z) (3) 

where 



J = 






Inxn 





Inxn 














Ofexfe 





















P = Mq-^\igi(q) 

i=l 

H(z) - \ (p + ]T XMq)) M- 1 (p + J2 X ^(l) ) + V(q) 



and, furthermore, the Euler- Lagrange equations for are equivalent to the 
generalized Hamiltonian system |#|). Eq. forms a constrained Hamiltonian 
system of the type |7p with constraint submanifold C a graph over (q,p), i.e., 
C := {(q,p, A) : A = X(q,p)} and restricted symplectic form i*cu = dq A dp. 

2 We do not use the term nonholonomic which is reserved for constrained mechanical sys- 
tems satisfying the Lagrange-d'Alembert principle, whose flow is not in general symplectic. 



2 



Proof. Introducing Lagrange multipliers Ai, . . . , Xk, the Euler-Lagrange equa- 
tions for Q are 

j t (V 4 F) V q F = 0, (4) 
9i{q)-4 = Q, i = i,...,k, (5) 

where 

1 k 
F(q,p, A) = -q T Mq - V(q) - £ \ l9l {q) ■ q. 

i=l 

Expanding out equation dil gives the Euler-Lagrange equations 



d_ 
db 



d 
db 



Mq-J2^9i(q)\ ~V 9 F = 

V i=l ) 

Mq - X ^(i) ) + \ W(?) + J2 ^ D 9Mq) = o (6) 



Define the conjugate momentum p € K™ using the standard Legendre trans- 
form 

k 

p~V 4 F = Mq-J2^9i(Q) (?) 

i=l 

so that 

<7 = M~ x ^p + ^2 \9i(q)J ■ (8) 

Using equations ^ and |8| in equation Q gives the expression for p 

k / k \ 

p = -VV(q) - ^ XiDg^M- 1 lp + J2 W«) J ■ (9) 

Define the Hamiltonian, again in the standard way, as H(q,p, A) := q ■ p — 
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F(q,q, A); explicitly, 
H = q-p-F 

k 

\q T Mq + V{q) + Y^X igi {q) 



q-p r 



q-\p- \Mq + AiSi(g) J + ^(?) 



2 

i=l 



9 • f Mg - ^ A l3l ( g ) - l -Mq + ^ Aifli(?)] + V(q) 

\ i=l ~ z=l / 



J ( p + E J Arl ( p + E «g) J + 



A calculation shows the equivalence of the right hand side of ^ and V p H; of 
the right hand side of Q and —V q H(q,p,X); and of constraints gi(q) ■ q = 
and = V x H(q,p, A). 

The constraints = V \H(q,p, A) are the following set of equations linear in 

A, 

(g x ■ M^gx ■■■ g x - M- V 





\g k ■ M~ x gi ■■■ g k - M~ x g k/ 

which has a unique solution for A for all q, p because the matrix is GM~ 1 G T 
where G is the k x n matrix whose «th row is gj . The assumption that the gi 
are linearly independent means that G has full rank k and hence that GM~ 1 G T 
is nonsingular. The constraints therefore have a unique solution for A that we 
write as A = X(q,p), that is, the constraint submanifold is a graph over (q,p). 
Differentiating these constraints with respect to t then yields ODEs for A, that 
is, the system (|3| has (differentiation) index 1. The symplectic form on C is 
\dz A Jdz = dq A dp. □ 



2 Symplectic integrators for generalized Hamil- 
tonian systems 

The Hamiltonian form (|3| suggests considering generalized Hamiltonian systems 
of the form Jz = VH(z), z £ K m , where J is a constant antisymmetric matrix, 
and we do not specify the constraints. Note that many kinds of constrained 
Hamiltonian systems (including those with holonomic constraints) can be writ- 
ten in this form; the constraint manifold C is constructed as the subset of initial 
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conditions for which the equations have a solution. In general, these equations 
may not have solutions for all initial conditions; in the extreme case J = 0, the 
equations are purely algebraic. However, it is easily seen that any solutions that 
do exist do preserve the (possibly degenerate) bilinear form u T Jv, for 

— u Jv — u Jv + u Jv = u H zz (z)v — v H zz (z)u = 

where H zz is the Hessian of H — this does not require invertibility of J. 

In the particular case of Proposition [T] the generalized Hamiltonian system 
that is obtained is equivalent to a canonical Hamiltonian system obtaine d by 
eliminating the Lagrange multipliers A. Let A = X(q,p) be the solution to (10 1. 
Then Hamilton's equations for H(q,p) := H(q,p, X(q,p)) are 

• dfl ( \ 
Qi = -t;— (q,P) 

OPt 

k 



OPi 

^-(q,p,Mq,p)) 

OPi 



-(q,P,Mq,p)) +52 ^r- A (q,p) 
dpi dX 3 dpi 



OH 

Pi = —k— \,q,p) 

dqi 



OH ~ x— -\ OH - OX 

-— (q,p,x{q,p)) -J2 {qipMq,p))-Q^{q,p) 
-^(q,p,Mq,p)) 



which, together with ^(q,p,X(q,p)) — 0, are equivalent to (|3j). That is, the 
two operations of eliminating the Lagrange multipliers and mapping the Hamil- 
tonian to its Hamiltonian vector field commute; this can also be seen abstractly 
by considering the symplectic manifold C with canonical coordinates (q,p), sym- 
plectic form dq A dp, and Hamiltonian i* H . 

The midpoint rule is known to be symplectic when the structure matrix J is 
invertible |3_. However, as for the continuous time case, J need not be invertible. 

Proposition 2. Any solutions of the midpoint rule applied to Jz = Vif pre- 
serve the 2-form \dz A Jdz. 

Proof. It must be shown that dz\ A Jdz\ — dzQ A Jdz$, where z\ is the result of 
applying the midpoint rule to z . The midpoint rule is 

J ^ = Vtf(^):=V^) (11) 
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where At is the time step and the method takes a point zq to z\ . 
Taking exterior derivatives of equation (111 yields 

Jdz\ — Jdzo + AtH zz (z) 

where H zz is the Hessian of H. Then 

J(dzi -dzo) = -AtH zz {z){dz Q + dzi) 

=> (dzo + dz\) A J[dz\ — dzo) — (dzo + dz%) A -AtH zz (z)(dzo + dzx) 

=> (dzo + dzi) A J(dz\ — dzo) = 
=>■ dz\ A Jdzi — dzf) A Jdz = 



This is an instance of the following more general result. 



□ 



Proposition 3. Any solutions of any symplectic Runge-Kutta method applied 
to Jz = V-ff preserve the 2-form ^dz A Jdz, where J is any constant antisym- 
metric matrix. 



Proof. The s stage symplectic Runge-Kutta method is 

s 

JZi = Jzq + At a,ij JFj 



4=1 



J Zl = Jz + At ^2 b j JF j i 

4=1 



where 



The coefficients of a symplectic Runge-Kutta method obey 

0. 



Taking the exterior product of equations (12), (13 1, and (14) gives 



Jdzg — JdZi — At fly JdFj 
4=1 

S 

Jdz\ = Jdzo — At frj Jdi^j 
JdFj — H zz dZj 



(12) 
(13) 

(14) 
(15) 

(16) 

(17) 
(18) 



From equation ( 18 I, 



dZj A JdFj = dZj A H zz dZj = 



(19) 



G 



Then 



dz\ A Jdzi — dz A Jdzo 

s 

= dz 1 A J(dz + At'^2b J dF j ) - dz A Jdz (using (l7|) 



3=1 



= — Jdz\ A dzQ + At bj dFj — dzo A Jdzo 



3=1 



= - Jdz + At ^2 b j JdF j A [ dz + At^2 fydFj - dz A Jdz a 



3=1 



3 = 1 



AtdZQ A J ^ b 3 dF 3 + At /2 b 3 dF 3 A JdZa + ^ h 3 dF 3 A J z2 b 3 dF 3 



3 = 1 



3 = 1 



3 = 1 



3 = 1 



-Aty^bjJdzo A dFj + At^bidFi A Jdz + At 2 ^ bjdFj A J^bjdFj 



3 = 1 



( = 1 



3 = 1 



3 = 1 



-At^bjJ ( <2Zj - At^Oj-idFj ) A dFj (using ( fl6| ) 



At ^b.dF.AJ dZ, - At ^ ctydFj (using ( fl6} ) 



= -Ai 2 



bjQjidFi A •/'// ', + V] h,(i, ,ili, A •/<// , 



(using (1191) 



+ A/ J ^ /).(//•', • ./^J /).(//•'. 

3=1 3=1 



= At 2 



^ - &j<iji - &iOij) dF* A JdFj 

1,3 = 1 



(using (15)) 



□ 



A full study of the geometry of the relations (zq,Zi) generated in Propo- 
sition [3] remains to be undertaken]^] Unfortunately, the relations (zn,zi) in 

3 The relations generated in Proposition [3] are a generalization of the Viterbo generating 
functions used in symplectic topology [JJ. These take the form S: QxR fc — > R; the submanifold 
p = S q (q,X), = S\(q,X) is Lagrangian in T*Q. The parameters A allow the representation 
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Proposition [3] do not yield good integrators for arbitrary J and H. For exam- 
ple, holonomic constraints can be specified as generalized Hamiltonian systems 
with H — H(q,p) + Xa=i Khi(q). In this case the midpoint rule, say, generates 
maps from all (qo,Po) to (qi,Pi) with the constraints satisfied at the midpoint. 
Not only is the phase space 'wrong', this method is known to be not convergent 
in general [2J. The situation is much better for the index one constraints of 
Proposition [T] 

Proposition 4. For the index 1 constrained problem of Proposition^ Propo- 
sition [5] yields integrators that are well-defined for sufficiently small At, con- 
vergent of the same order as the Runge-Kutta method, preserve the constraint 
submanifold, and preserve the symplectic form on the constraint submanifold. 

Proof. In this case the constraint part of the Runge-Kutta equations read 

= ^ x H(Q i ,P i ,A i ), i = l,...,k, 

A i— ^- = y £b t V x H(Q i ,P i ,A i ) 



Therefore the Lagrange multipliers A* at each stage are given by the exact La- 
grange multipliers evaluated at (Qi,Pi), i.e. Aj ~ X(Q i ,P i ), and Ai is arbitrary. 
For convenience, we add the extra equations Ao = X(qo,po), Ai = X(qi,pi) which 
do not affect the method at all. The resulting method is equivalent to that ob- 
tained by eliminating the Lagrange multipliers in the Hamiltonian, applying a 
symplectic Runge-Kutta method, and lifting back to the constraint manifold by 
A = X(q,p). It is therefore well defined for sufficiently small At and convergent 
of the same order as the Runge-Kutta method. Because ^dz A Jdz — dq A dp, 
the symplectic form dq A dp is preserved on the constraint manifold. □ 



3 General constraints 

Under certain conditions, namely that the Legendre transform that defines the 
conjugate momenta must be invertible to give q, Proposition [T] can be gener- 
alized to allow a general Lagrangian and general constraints. A very thorough 
geometric treatment of this type of constraint, applying the Gotay-Nestor ge- 
ometric version of the Dirac-Bergmann constraint algorithm, can be found in 

Proposition 5. If the Legendre transform mapping (q,q,X) —> (p,q,X) given 



in equation (21) is invertible then the Euler-Lagrange equations for the action 



S{q)= / L{t,q,q)dt 



of larger classes of Lagrangian submanfolds than the standard generating function S(q) which 
generates p = S q (q) which is necessarily a graph over Q. 
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subject to the constraints gi(q, q) = 0,i = 1, . . . , k are equivalent to the general- 
ized Hamiltonian system 

Jz = VH(z) 



(20) 



whe 



"■ri x n 

o 

p 





feX fc 



(21) 



p = VgF(q,q,X) 
H(z) = q-p-F(q,q,X) 
k 

F(q,q,X) = L(t,q,q) - 2J Xi9i (q, q) 



(Here the Lagrange multipliers in the (q,q) formulation are determined by ( pj 
^25^ below.) If, in addition, the matrix G(q,q) given by Gij = dgi(q,q)/dqj has 
full rank k for all q, q then the system of Eq. J[2(ty has index one, i.e., can be 
solved for A = X{q,p). 



(22) 



Proof. The Euler-Lagrange equations are 

d 



dt 



(V 4 F) - V q F = 0, 



gi(q,q)=0, i = l 



, . , . , k, 



where 



F(q, q, A) = L(t, q,q)-^2 A *ft(<7, <?)■ 



Define the conjugate momentum p as 



By assumption equation ( 26 ) can be rearranged to give q 

q = f{q,pA)- 



Using equation ( 26 ) in equation ( 23 I gives the expression for p 



p = V 9 F. 



Define the Hamiltonian as 



H(q,p,\) :=q-p- F(q,q,X) 

= P- f(q,PA) - F(q,q,X) 



using (27 1 



(23) 
(24) 

(25) 

(26) 

(27) 

(28) 

(29) 
(30) 
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Then (where f p is the Jacobian derivative dfi(q,p,X)/dpj, etc.) we have 



\7 p H 



V p P 



f + fi 
= f + f P P ~ / P V 4 P 
= f + f P P ~ f P P 

= <i 

V q H = f q p-f g V q F-V q F 

= fqP - fqP - V g F 
= ~P 

V A ff = fxp - fx^ q F + g 
= f\P - fxP + 9 
= 



using (26) 



using (27) 



using (26) 



using (28) 



using (26) 



using (24) 



(31) 



(32) 



(33) 



establishing the equivalence of the Euler-Lagrange equations to ( 20 1 . The con- 
straints in the Hamiltonian formulation are 

= VH(q,p, X)=g(qJ(q,p, A)) 

and their Jacobian derivative with respect to A is the matrix 

% dqj 



E 



dqj dX k 



The first factor is G. The second factor is the derivative of the inverse Legendre 
transform (assumed invertible) with respect to A. The forward Legendre trans- 
form is p = + J2i=i ^i^ii anc ^ ^ s derivative with respect to A is G T . By the 
chain rule, if G has rank k then the Jacobian is nonsingular for all q, q and the 
constraints have a unique solution A = X(q,p) for all q,p, that is, the system 
has index one. □ 

The integrators we considered in Section [2] work for any index one system. 

Proposition 6. For any constant antisymmetric J, if the generalized Hamil- 
tonian system Jz = ViJ(z) can be solved for X — X(q,p) where z — {q,p 1 A) are 
Darboux coordinates for J, then any symplectic Runge-Kutta method applied 
to this systems yields constrained symplectic integrators convergent with their 
classical order. 

Note that the assumptions are satisfied if |-Haa| 7^ 0. The constraints may 
be nonlinear in A, and need not be solved analytically; the entire Runge-Kutta 
system for (Qj,Pj,Aj) can be numerically solved simultaneously. 

Proposition[5]can be generalized further, to any singular Lagrangian L(q, q, A), 
and still further to Lagrangians L(z,z) where \La\ = 0, but the required non- 
degeneracy assumptions are not as geometrically transparent as those in Propo- 
sition |5l 
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4 Sub-Riemannian with holonomic constraints 



Proposition [T] allowed a variational problem with subRicmannian constraints 
to be converted into an unconstrained Hamiltonian system that can be solved 
using the symplcctic midpoint rule. The proposition in this section shows that 
if holonomic constraints are added to the original variational problem, then 
the resulting Hamiltonian system is a simple holonomically constrained system. 
This system can be solved by a symplectic method such as rattle [5]. 

Proposition 7. Let M be a symmetric nonsingular nxn mass matrix, V : 
W l — >• R a smooth potential, gi : R™ — > M. n , i = 1, . . . , k be k smooth functions , 
and q be a smooth extremal with fixed endpoints for the functional 

S(q) = jf 1 L(t, q, q)dt = jf 1 (^q T Mq - V(q)^j dt (34) 

subject to the velocity constraints gi(q) ■ q = 0, i = 1, . . . , A; and the holonomic 
constraints hi(q) — O.i — 1, . . . , I. Then 

Jz = VH{z) (35) 

where 



J = 



( 














(A 


In x n 















p 








Ofexfe 







A 


V o 








Oix 









k 

p = Mq-Y x i9i{q) 

i=l 

H(z) = \(p + Y X ^)) M_1 (p + E X ^(l)) + V(q) + £ X^hiiq) 

\ i=l / \ i=l / i=l 



and, furthermore, the Euler-Lagrange equations for ( 34 ) are equivalent to the 
generalized Hamiltonian system (35). If, in addition, the velocity constraints 



are linearly independent for all q, then Eq. (35) is equivalent to a canonical 
holonomically constrained Hamiltonian system. 



Proof. As in Proposition [T] the extended Lagrangian F, the conjugate momenta 
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p, and the Hamiltonian H(q,p, A, X h 



) are defined by 



F 




P 



VqF = Mq-J2*i9i(q), 



i=l 



H 



q-p-F. 



The rest of the proof is a calculation along the same lines as for Proposition 



Proposition 8. Subject to standard nondegeneracy assumptions on the Hamil- 
tonian, the following algorithm yields a convergent, second order integrator that 
is symplectic on the constraint manifold defined by the (primary) holonomic con- 
straints and the secondary constraints induced by them: (i) apply RATTLE using 
the holonomic constraints; (ii) in the inner step of RATTLE, when a time step of 
the unconstrained system is required, apply the midpoint rule to the generalized 
Hamiltonian system with Hamiltonian H(q,p, A, 0). 

Proof. Eliminating the velocity constraints by solving for the Lagrange multi- 
pliers yields a standard holonomically constrained system. Applying RATTLE 
(with the midpoint rule in the inner step) to this system yields a convergent 
second order integrator on the constraint surface. Applying the midpoint rule 
in the inner step is equivalent to applying the midpoint rule to the generalized 
Hamiltonian system with Hamiltonian H(q,p, A, 0). □ 

5 Example: Sub-Riemannian geodesies 

The motion of a two-wheeled vehicle with a front steering wheel and a non- 
steering back wheel, moving on a smooth surface, will be modelled. We consider 
the two wheeled vehicle shown in Fig. [I] with length L, back wheel at (z,w), 
and front wheel at (x,y). The front wheel is at an angle <fi and the vehicle is at 
an angle 9. 

If the speed of the front wheel is v, its velocity of the front wheel must obey 



D 



□ 



x = v cos <fi 
y = v sin <p 



Eliminating v, the velocity of the front wheel obeys the constraint 



x sin </> — y cos <f> = 



(36) 



Similarly, the velocity of the back wheel obeys the constraint 



i sin 9 — w cos 9 = 



(37) 



12 




(z,w). 



Figure 1: A two wheeled vehicle with the front wheel at an angle 
entire vehicle on an angle of 9 



and the 



We can eliminate equation (37 1, and thus the variables z and w, using the 
distance between the two wheels which relates the four variables. Notice that 



x — z = L cos 9 
y — w = L sin 9 

which after taking time derivatives gives 

i = x + L9 sin 9 
w = y — L9 cos 9 



which substituted into equation ( 37 1 gives 



x sin 9 — y cos 9 + L9 = 



(38) 



The constraints given in equations ( 36 1 and ( 38 1 can be written as 



9i{l) ■ Q = 0,i = 1,2 



where q — (x, y, 9, <j)) T and 



9i (<?) = (sin (j), - cos <f), 0, 0) T 
Qi (<z) = (sin 9 : — cos 9,L,0) T 



We take the Lagrangian to be 
1 



L = 



+ y 2 + a9 2 +p<p 2 \-V(x,y) 



(39) 
(40) 



(41) 



where the potential V(x, y) is the (scaled) height of the surface, giving an index 
one system as in Proposition [l] That is, we are calculating geodesies (in the 
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Figure 2: Snapshots at every 600 step of the bicycle starting at 
{x,y,6,<j>,p x ,Py,peiP4>) — (0.15, 0,7r/4,0, 1,0, 0), with At = 0.001, and trav- 
elling in a potential V(q) = 1 + cos(r). At about t — 7.8 and t = 15.6 the 
steering wheel is aligned with the bicycle and the bicycle changes direction and 
retreats rapidly. 



case V = 0) of the subRiemannian metric defined by Eqs. (41 1, (36) and (38 1. 
We use the midpoint rule. 

The first four tests use the potential V(q) — — cos r. where r is the midpoint 
of the vehicle. The first test is to compute a simple trajectory of a bicycle of 
length 0.3, starting with the centroid of the bicycle at the origin. The bicycle 
is given a small push by setting the initial generalized momenta p y ^ 0. The 
initial conditions are {x,y,6,(j),p Xi p y ,pe,p^) = (0.15, 0, 7r/4, 0, 1, 0, 0). Fig. [2] 
shows the bicycle at regular snapshots through time. The background colours 
show the potential: red is high, and blue is low. Note that that at about t = 7.8 
and t — 15.6 the steering wheel is aligned with the bicycle and the bicycle 
changes direction and retreats rapidly. 

The second test is to check the order of the method by plotting, in Fig. [3] the 
error of various solutions. A reference trajectory with At = 10~ 4 and final time 
10s is computed, and trajectories with bigger time steps are compared to it. The 
slope of the error line shows that the method is second order. A comparison 
to a highly accurate reference solution calculated with matlab's odel5s and 
tolerance 10 -12 also shows the second order accuracy of our method (Fig. El. 

Our third test is to numerically test the symplectic condition by evaluating 
the symplectic bilinear form for a number of steps. See Fig. [5j The reference 
trajectory is the one used in the error plot with with At = 0.1. Two nearby 
trajectories are also calculated to allow an estimate of the change in dqAdp. The 
calculated change is near roundoff indicating that the integrator is symplectic. 

The fourth test is to plot, in Fig. [6j the energy error for At = 0.01 and 
At = 0.001. From the figure the energy errors appear to be bounded, as expected 
for a symplectic integrator. 

The free motion case (V(q) = 0) has two simple solutions that are relative 
equilibria for the translation and rotation symmetries of the problem, namely 
straight line and circular motion. The first trajectory to check is a straight line. 
If the bicycle starts with 9 = <fi = 0, and there is no potential field, then bicycle 
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Er ror compared to At = 10 4 solutio n 




icr 4 io" 3 10" 2 10 1 

At 



Figure 3: The error of three runs compared to a reference solution with 
At = 1 x 1CP 4 . The initial conditions were (x,y,0,4>,Px,Py^P8,P(f>) — 
(0.15,0,0, f,0, 1,0,0). The method is order At 2 



, Error compared to odel5s 

io -3 

10* 

o 

L1J io 6 

IO" 7 
IO" 8 

io 4 io 3 icr 2 io 1 

At 

Figure 4: Error compared to MATLAB ode 15s solution. Initial conditions were 
[x, y,9, 4>,p Xl p y ,pe,P4,) = (0.3, 0, 0, n, 0, 1.09, 1), with A = (0, -0, 3), final time 
10s, and travelling in no potential field. The solution agrees with the MATLAB 
solution. 
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, ie-i7 Bilinear form 



0.5 




-1.0 



1 S.O 0.2 0.4 0.6 0.8 1.0 
Step le2 

Figure 5: The symplectic bilinear form error evaluated for 100 steps with At = 
0.1. The trajectories cj> Zi (ji) for three nearby initial conditions, zn, Zx, and z%, 
were calculated. The change in the symplectic form was estimated as u^Jv n — 
UqJvq where u n — <f> Zl (n) — <fr z „(n) and v n — <j> Z2 (n) — <j> Zo (ri). The initial 
points were z = (0.15, 0, 0, f , 0, 1, 0, 0), z x = (0.1500000085,1 x 10- 8 ,-l x 
10~ 8 , 0.785398145543,-1 x 10" 8 , 0.99999998, 1 x 10~ 8 ,-1 x 10" 8 ) and z 2 = 
(0.1500000115,-1 x 10- 8 , 1 x 10~ 8 , 0.785398181251, 1 x 10~ 8 , 0.99999998,-1 x 
10~ 8 ,-1 x 10~ 8 ). 



Energy error 




Time 



Figure 6: The energy error over time. This is the energy at each step minus 
the initial energy. The bicycle is trapped in the potential bowl — V(q), and the 
energy error does not show a linear growth in time. 
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Trajectory snapshot every 600th step 
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Figure 7: Snapshots at every 600 step of the bicycle starting at 
(x,y,9, <p,p x ,p y ,pe,Pci>) — (0.15,0,0,1,0,0,0), with At — 0.01, and travelling 
in no potential field. The bicycle stays straight. The trajectory is unstable and 
eventually wanders from a straight path. 



should remain travelling in a straight line. 

Let 6 = <f) = 0, x = 1, and y = 0. The constraints in equations (361 and (381 
are satisfied. Equation (|8| gives the initial generalized momenta values: all are 
zero except p x = 1. In Fig. [7] this simple trajectory of the bicycle is confirmed. 
For the circle, let 9 — at, <p = at + | , x = — csin(0), and y — ccos(0). There 



are two constants, a and c, to be determined. Equation (36 1 gives A = (1, — c), 
and equation (38 1 gives aL — c. Using these values in equation d8| gives the 
initial generalized momenta values: (p x ,Py,Pe,P<f>) = (0,0, a(l +Lr),a). For 
this simple trajectory a is chosen to be 1. In Fig. [8] the circle trajectory of 
the bicycle is confirmed. If the trajectory is computed for larger times the 
bicycle leaves the circle; the solution appears to be unstable, but, interestingly, 
appears to repeatedly return to the circular orbit, indicating a possible relative 
homoclinic structure of this problem. This is also suggested from the evolution 
of the Lagrange multipliers, A 



(1, — c) being shown in Fig. 
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6 Example: the Heisenberg problem 

A previous study of geometric integrators for subRiemannian variational prob- 
lems used a discrete variational approach to obtain constrained symplectic inte- 
grators pp. Our approach, applying symplectic integrators to the Hamiltonian 
formulation, yields geometric integrators with the same geometric properties, 
but uses standard integrators that allow any order with standard implementa- 
tions, and does not require an approximation of q, that is, it naturally yields 
first-order trajectories in (q,p) instead of second-order trajectories in q. 

We repeat the numerical illustration of fl, pg. 12], the Heisenberg problem, 
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-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 
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Figure 8: Snapshots every 5 step of the bicycle starting at 
(Xjy^O^rp^PxjPyjPejPc/,) = (0.3, 0, 0, %, 0, 1.09, 1), with At — 0.1, and trav- 
elling in no potential field. The bicycle stays in a circle for many revolutions 
(not shown for clarity), but the trajectory is not stable, so eventually wanders. 



Circle, At =o.oi 
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Figure 9: Some of the phase space variables for the bicycle trajectory starting 
in a circle. The solution for pg and p^, suggest a relative homoclinic orbit. 
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Figure 10: Some of the phase space variables for the bicycle trajectory starting 
in a circle. These angle variables are monotonically increasing. 
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Lambda's over time 



- K 




"S.O 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 
Step le3 



Figure 11: The error in the A's for the circle trajectory. Note that there are 
many more steps than shown in Fig. [8] This shows that the A's periodically 
return to their values for a circle. 



using our approach. This is to find the extremal q(t) — (x(t), y(t), z(t)) of 



S(q) 



L(t,q,q)dt = 



q 1 q - V(q) dt 



subject to the constraint g(q) ■ q = 0, where g{q) = (—y,x, 1). 
Equation ^ gives q 




Using equation (|9| the p can be written 
= -W(q) - A 





and we have the constraint g ■ (p + Xg) = 0, which gives 

g-p 

9-9 



X 



(42) 



(43) 



A simple trajectory starting with the same initial conditions as in [TJ 
pg. 15] is shown in Fig. 12 Their initial conditions are (x, y, z, x, y, z, X) — 
(0,0,0,0.1,0.3,0,1), which when converted to generalized momenta variables 
are (x,y, z,p x ,p y ,p z , X) — (0,0,0,0.1,0.3,1,1). Qualitatively the results look 
like [TJ pg.14]. 



7 Discussion 

We have constructed symplectic integrators for a different class of constrained 
Hamiltonian systems than the holonomic constraints most commonly consid- 
ered in the literature. The class includes important practical problems arising 
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Heisenberg system 




Figure 12: The Heisenberg example starting at (x,y, z, x,y, z) — 
(0,0,0,0.1,0.3,0) or (x,y, z,p x ,p y ,p z , X) = (0,0,0,0.1,0.3,1,1). Qualitatively 
the results look like [TJ pg.14]. 



in subRieniannian geometry. We have restricted our attention to symplectic 
Runge-Kutta methods; a generalization to partitioned methods in which differ- 
ent Runge-Kutta coefficients are used for q, for p, and for A is straightforward. 
In other work [5] , we reinterpret these methods as an instance of rattle in an 
extended phase space; that point of view also suggests different generalisations. 
However, the nondegeneracy conditions are essential for the method, indeed, 
for the entire approach, to work. It is not clear to what extent the approach 
can be extended to handle more general constraints, for example, to the system 
Ji = Vff + AV5, where the constraint submanifold g(z) — is symplectic. 
No symplectic, constraint-preserving method is known for this problem. As 
remarked before Proposition [4] a full study of the geometry of the relations 
(z , Z\) generated in Proposition [3] remains to be undertaken. Any solutions are 
symplectic, so this gives access to a much larger class of symplectic maps than 
do traditional generating functions. Note that new variables (analogous to A) 
can be added as needed to generate larger classes of maps. 
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